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P.  O.  BOX  8172  ALBUQUERQUE,  NEW  MEXICO  87198  (505)W¥4JW  262 -0440 

15  October  1982 

Major  Carl  Edward  Oliver 
AFOSR 

Directorate  of  Mathematical  and  Information  Sciences 

Building  410 

Bolling  AFB,  DC  20332 

SUBJECT:  FIRST  RESEARCH  PROGRESS  AND  FORECAST  REPORT  (10/15/82) 

FOR  CONTRACT  F49620-82-C-0064 ,  "ADAPTIVE  GRID  GENERATION 
USING  ELLIPTIC  GENERATING  EQUATIONS  WITH  PRECISE 
COORDINATE  CONTROLS" 

Dear  Major  Oliver: 

This  letter  constitutes  the  required  report  on  the  subject 
contract,  for  the  5-month  period  from  the  inception  of  the 
contract  on  15  May  1982. 

Generally,  I  am  extremely  pleased  with  the  progress 
accomplished  and  the  forecast,  detailed  in  the  following 
topical  areas. 

(1)  2D  GRID  GENERATION 

The  codes  for  the  generation  of  the  2D  grid  equations , 
the  solution  of  the  finite  difference  equations ,  and  the 
verification  process  have  all  been  successfully  developed. 

See  #  2  below  for  details. 

(2)  3D  GRID  GENERATION 

We  have  completed  development  of  codes  for  .the  transformation 
of  general  second-order  equations  in  3D  into  general  nonorthogonal 
coordinates,  for  the  solution  of  these  equations  by  hopscotch 
SOR,  for  rigorously  verifying  the  algebra  and  coding  involved, 
for  establishing  the  grid  generation  equations  and  solving 
them,  and  for  rigorously  verifying  these. 

Using  Symbolic  Manipulation* on  a  VAX  computer,  code  was 
written  to  transform  the  general  second-order  PDE  to  general 
nonorthogonal  coordinates.  The  result  is  a  19-point  operator, 

-  83  .!  ?  o94 

*The  Symbolic  Manipulation  work  was  done -under  related  contract 
from  the  U.S.  Army  Research  Office,  and  involved  Dr.  S.  Steinberg. 
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since  the  8  corner  points  on  the  cube  of  the  3D  operator  are 
identically  zero.  A  Fortran  subroutine  was  written  to  generate 
the  3D  arrays  for  the  finite  difference  coefficients .  These 
were  solved  by  hopscotch  SOR  code. 

The  validation  procedure  for  the  "hosted  equations"  (in 
this  case,  the  variable-conductivity  electric  field  equation) 
consisted  of  testing  the  truncation-error  convergence  of  the 
solution.  An  inverse  procedure  was  devised,  in  which  the 
continuum  solution  was  specified,  chosen  so  as  to  possess 
enough  structure  to  exercise  all  the  derivatives  of  the  operator 
and  all  finite-difference  truncation  errors.  For  the  second- 


order  operator  and  second-order  accurate  finite  difference  forms, 

the  solution  was  specified  as  sol  =  x3y4z5.  The  transformation 

used  involved  the  hyperbolic  tangent  of  the  product  of  all 

three  transformed  coordinates,  =  x^  +  tanh(d^xyz) ,  where 

x^  =  x,  X2  =  y,  Xg  =  z.  The  hosted  equation  in  the  original 

coordinates  was  L($)  s  v-aV(j>  =  q.  The  variable  conductivity 

0  was  also  chosen  to  give  significant  structure  to  the  matrix, 

with  o=o  sin(xyz/x  y  z  ) .  The  non-homogeneous  part 
o  max  max  max 

q  was  chosen  so  as  to  give  the  desired  solution,  i.e.  q  =  L(sol). 
This  highly  structured  problem  was  then  fed  to  the  Symbolic 
Manipulation  code,  and  the  matrix  problem  generated  by  it  was 
solved  numerically.  (In  2D,  we  used  the  GEM  spatial  marching 
codes,  and  in  3D,  we  used  the  hopscotch  SOR  solver.)  By 
monitoring  the  truncation  error  as  the  grid  was  refined  from 
53,  93,  173,  333,  we  verified  the  transformation,  the  finite 
difference  forms  (validating  the  0(A2)  accuracy)  and  the  itera¬ 
tive  solution  procedure.  As  expected  theoretically,  the  value 
of  C  =  A2  *TE,  where  TE  is  the  maximum  truncation  error  in  the 
mesh,  becomes  constant  as  the  mesh  is  refined.  The  size  of 
C  depends  on  the  grid  stretching  parameters  d^  and  oQ,  but  the 
entire  method  remains  0(A2)  accurate. 

A  similar  procedure  was  followed  for  the  grid  generation 
equations  themselves,  although  this  was  more  complicated. 


A  "solution"  for  the  grid  transformation  problem  L(£^)  ■ 

was  chosen  as  £.  -  x .  +  tanh(D.xyz) ,  where  the  D.,  are  grid 
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stretching  parameters,  not  the  same  as  the  ch  used  earlier. 

(Note  the  roles  of  and  5^  are  reversed  from  the  earlier 
case  of  the  hosted  equation  validation.)  The  non-homogeneous 
terms  are  then  chosen  to  give  the  desired  solution, 

=  L(g^)  where  L  is  the  Laplacian  in  x,y,z.  With  specified 
at  each  mesh  point  in  the  transformed  domain,  the  x.  had  to 
be  obtained  by  solution  of  a  3x3  transcendental  equation  system. 
Because  of  round-off  sensitivity  (see  below)  this  was  accomplished 
with  a  3x3  Newton-Raphson  procedure,  converged  to  essentially 
single  precision  on  the  computer.  The  numerical  solution 
procedure  for  the  grid  generation  equations  was  also  more  involved 
then  the  hosted  equation  solution,  since  it  involved  a  nonlinear 
system  of  three  coupled  equations  for  x,y,z.  This  solution  was 
obtained  by  outer  (nonlinear)  Picard  iterations  and  inner  (linear) 
hopscotch  SOR  solutions.  The  result  and  interpretation  is 
the  same  as  the  hosted  equation  situation,  with  constancy  of 
C  =  a2-TE  validating  the  transformation,  the  finite  difference 
forms,  and  the  solution  procedure,  with  the  grid  varied  from 
53  ,  93  ,  173 ,  333 . 

The  calculations  were  performed  on  a  32 -bit  computer,  a 
VAX  780,  with  approximately  8  significant  decimal  figures  of 
accuracy.  The  round-off  problem  only  slightly  obscured  the 
truncation-error  convergence  testing  for  the  hosted  equations , 
but  was  a  serious  problem  for  the  grid  generation  equations . 

A  false  indication  of  a  persistent  error,  which  stopped  progress 
for  the  better  part  of  a  month,  was  eventually  traced  to  the 
evaluation  of  the  Laplacian  of  the  solution,  which  involved 
three  summations  of  terms  containing  sech2  and  tanh.  We  are 
still  investigating  the  reason  for  the  error,  but  it  was  cured, 
and  the  validation  procedure  was  made  successful ,  with  a  minor 
grouping  of  terms  in  the  Laplacian  subprogram. 

Although  this  round-off  error  problem  would  have  been 
nonexistent  on  a  CDC  or  Cray  computer,  we  are  still  convinced  that 
the  excellent  interactive  editing  and  running  capabilities  on 
the  VAX  computer  was  essential  to  the  rapid  development  of  these 
codes.  In  fact,  the  Symbolic  'Manipulation  code  used  is  available 
only  on  the  VAX  computers.  However,  even  the  other  code  components, 
involving  more  traditinal  numerical  methods ,  were  developed 
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much  more  rapidly  than  they  would  have  been  on  a  mainframe 
computer.  We  estimate  a  factor  of  3  savings  in  code  development 
time,  with  better  structured  code  resulting  also. 

(3)  TRUNCATION  ERROR  VERIFICATION  r 

— < 

The  question  of  the  accuracy  of  grid  transformation  techniques 
has  been  confused  recently.  Mas tin  and  others  have  gone  through  : 

lengthy  analyses  and  shown  that  second-order  accuracy  is  lost 
for  severe  distortions.  These  analyses  are  in  disagreement  r4 

with  previously  published  results,  including  some  by  the 
present  P.I.  Mas tin  has  some  experiments  which  he  calls  on  ! 

to  validate  his  results .  The  numerical  experiments  show  that 
the  truncation  error  increases  as  the  grid  distortion  increases .  f 

This  work,  and  similar  analyses,  were  the  subject  of  private 
discussions  at  the  Nashville  conference  on  Grid  Generation  held 
in  April  1982.  I  am  of  the  opinion  that  these  analyses  are 
incorrect,  and  that  the  numerical  experiments  do  not  confirm  p* 

or  contradict  the  analysis .  Everyone  agrees  that  severe  distor¬ 
tion  will  increase  the  truncation  error,  but  the  question  of  : 

the  order  of  the  accuracy  can  only  be  settled  by  doing  the  kind  ■ 

of  systematic  truncation-error  testing  which  we  have  done  in  the  jp 

present  work.  The  present  results  show  clearly  that  second-order  i 

accuracy  is  maintained,  even  for  strong  grid  distortion. 

(4)  2D  ADAPTIVE  GRID  FOR  E-FIELD  CALCULATIONS 

A  procedure  and  code  has  been  successfully  developed  for  jr 

solution-adaptive  gridding  for  E-field  calculations  in  2D.  1 

The  calculations  are  being  performed,  under  separate  funding, 
for  the  Air  Force  Weapons  Laboratory  in  support  of  a  laser  j 

development  project.  The  solution-adaptive  procedure  was  developed  |H 

under  present  AFOSR  funding. 

This  E-field  calculation  differs  from  most  fluid  dynamics 
calculations  in  that  the  maximum  value  of  the  electric  field 
occurs  on  the  boundary.  (This  is  rigorously  true  for  the  linear 
case  with  constant  conductivity,  and  is  expected  in  all  practical 
cases.)  Thus,  the  grid  adaptivity  does  not  concern  the  grid 
point  location  at  interior  points ,  obtained  with  tailoring  of 
the  nonhomogeneous  terms  in  the  tranformation  equations,  but 
the  location  of  the  grid  points  on  the  boundary.  This  aspect 
of  the  grid  transformation  is  also  of  interest  to  fluid  dynamic 
calculations,  however,  and  in  fact  has  been  largely  neglected. 


□  □ 


In  previously  published  papers,  the  grid  distribution  along  the 
boundaries  has  been  obtained  in  an  ad  hoc  manner,  with  no  adaption 
used. 

In  the  presently  developed  procedure,  we  obtain  a  first  solu¬ 
tion  on  a  grid  with  the  boundary  distribution  determined  by 
the  local  curvature  of  the  surface.  (Default  method  consists  of 
equidistribution  of  points  over  arc  length,  but  this  can  be 
weighted  to  pack  more  points  near  local  surface  curvature.  With 
the  solution  adaptive  procedure,  we  find  that  equidistribution 
over  arc  length  is  adequate.)  The  second  grid  is  obtained  with 
boundary  points  packed  more  closely  near  the  maximum  E  values, 
using  an  empirically-determined  (user-specified)  weighting 
factor.  The  process  can  be  iterated  to  a  convergence  of  the 
grid,  but  in  practice  a  single  adaptive  step  is  sufficient  to 
attain  the  increase  in  resolution  desired. 

Overall  accuracy  (e.g.  in  :in  L2  norm  on  4>  or  E)  is  not  the 
objective  in  this  procedure,  but  rather  the  resolution  of  the 
maximum  value  of  E.  It  is  this  maximum  E  which  limits  the  power 
output  of  the  laser,  since  it  determines  arcing.  Although  the 
first  (non-adapted)  grid  may  give  E  values  which  are  quite 
accurate,  say  to  0.17o,  the  adaptive  step  increases  the  value 
of  the  maximum  calculated  E  value  by  as  much  as  20%  in  geometries 
considered  so  far.  It  appears  that  this  procedure  can  reduce 
the  computing  time  to  achieve  a  given  level  of  resolution  of 
maximum  E  by  a  factor  of  2  to  5,  compared  to  brute-force  grid 
refinement. 

Although  conceptually  simple,  the  implementation  of  this 
procedure  involves  a  subtle  procedure  of  what  we  have  called 
"projective  interpolation".  The  "true"  or  continuum  surface 
is  defined  by  a  set  of  discrete  coordinates  and  an  associated 
interpolation  rule;  we  refer  to  the  discrete  coordinate  and 
the  interpolation  rule  collectively  as  GQ.  The  initial  grid 
is  layed  down  by  some  rule,  say  the  equidistribution  in  arc 
length,  by  reference  to  G  ,  and  is  "on"  G  .  But  the  discrete  y.(x 
differ; the  domains  of  definition  are  not  the  same.  For  example, 

G  may  be  defined  by  1000  points  and  a  cubic  spline  interpolation 

w 

rule,  whereas  the  initial  grid  G^  might  consist  of  31  points 
on  the  boundary,  none  of  which,  other  than  the  end-points, 

need  correspond  to  any  of  the  1000  points  of  Gc-  When  the 
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second  (solution-adapted)  grid  G2  is  defined  by  resolution 
of  the  E-field,  which  has  been  calculated  not  on  Gc  but  on 
Gp  the  resulting  grid  points  will  be  "on"  G^  (say,  in  the 
sense  of  being  linear  interpolants  between  points  of  G^) 
but  not  generally  "on"  Gc-  (This  is  true  even  if  Gc  was 
defined  with  a  linear  interpolation  rule,  as  long  as  the 
defining  points  of  Gc  were  not  the  first  grid  G^.  This  is 
the  only  practical  case  to  consider.)  In  order  to  assure  that 
the  solution  adaptive  grid  is  actually  "on"  Gc,  a  combination 
of  interpolation  and  projection  is  necessary.  This  gets 
particularly  tricky  near  acute  angles  in  the  boundary,  which 
are  not  of  practical  interest  in  the  laser  electrodes  but 
could  be  of  interest  in  aerodynamics  problems,  where  the 
analogous  grid  adaption  would  be  motivated  by  the  desire  to 
accurately  resolve  separation/reattachment  points.  Even 
without  the  aci:te  angle  situation,  the  projective  interpolation 
problem  is  difficult.  The  details  of  our  "solution"  will 
not  be  given  here,  but  suffice  it  to  say  that  the  solution  is 
heuristic.  Note  that  when  the  interpolated  function  and  the 
defining  function  have  different  domains  of  definition,  we 
are  out  of  the  realm  of  analysis.  (For  example,  their  is  no 
mean  value  theorem.) 

We  have  given  some  consideration  to  the  analogous  3D  problem 
and  have  made  some  progress ,  but  have  not  completed  this  work 
nor  done  any  testing. 

(5)  3D  SURFACE  GRID  GENERATION  AND  INTERIOR  GRID  CONTROL 

Aside  from  the  solution-adaptive  problem,  and  the  associated 
projective  interpolation  problem  mentioned  above,  the  simple 
non- adaptive  generation  of  the  3D  grid  on  non-planar  surfaces 
is  a  difficult  problem.  Ruppert  of  Boeing  considers  this  the 
most  difficult  part  of  3D  aircraft  calculations.  P.D.  Thomas 
has  a  good  method  which  depends  on  projections  from  the  aircraft 
boundary  onto  some  outer  surface.  Although  successful  for  a 
wide  range  of  geometries,  it  will  fail  for  difficult  shapes  where 
the  projection  is  multivalued. 

We  have  made  some  progress  on  a  method  of  generating  the 
surface  grid  using  modified  3D  grid  generating  equations  solved 
in  the  non-planar  surface.  The  concept  would  accomplish 
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the  surface  grid  generation  in  a  smooth  but  otherwise  general 
surface;  e.g.  zsurface  could  be  a  multi-valued  function  of 

xsurf ace ’  ysurface'  We  have  not  advanced  to  the  coding  and  testing 
stage  as  yet. 

Likewise,  we  have  not  coded  and  tested  the  methods  for 
internal  grid  point  precise  control.  This  will  be  addressed 
shortly,  for  the  2D  problem  first. 

(6)  PUBLICATIONS 

We  will  present  an  invited  paper  entitled  "Symbolic  Manipulation 
and  Computational  Fluid  Dynamics"  at  the  AIAA  Computational 
Fluid  Dynamics  Conference  to  be  held  July  1983  in  the  Boston 
area.  The  primary  funding  for  this  work  was  provided  by  the 
Army  Research  Office,  but  the  paper  will  include  much  of  the 
work  done  under  the  present  AFOSR  contract  as  well. 

fbr  the  same  meeting,  I  will  be  submitting  a  contributed 
paper  on  the  adaptive  surface  gridding  problem  described  above. 

If  progress  is  rapid,  I  may  also  submit  another  paper  on  the 
3D  surface  grid  generation  problem.  Also,  I  will  submit  an 
Open  Forum  paper  on  the  truncation-error  validation  problem, 
and  try  to  lay  to  rest  this  question  of  the  order  of  accuracy 
of  highly  distorted  grids.  This  result  would  then  be  submitted 
to  the  AIAA  JOURNAL  as  a  Technical  Note. 

I  have  been  invited  to  serve  as  the  Editor  for  a  Special 
Volume  of  COMPUTER  METHODS  IN  APPLIED  MECHANICS  AND  ENGINEERING 
dedicated  to  Grid  Generation.  This  issue  will  probably  appear 
in  early  1984. 

(7)  BUDGET 

All  expenditures,  both  direct  charge  and  overhead  items, 
are  within  the  budget  for  this  contract.  We  do  not  anticipate 
any  changes  in  the  expenditures  or  time  allotments  from  the 
contract  budget  and  SOW. 

If  further  details  are  needed  on  any  of  these  topics, 
please  do  not  hesitate  to  contact  me  at  the  above  address  or 
telephone  number. 


Respectfully, 


Dr.  Patrick  J.  Roache,  Prin.  Inves 


